arXiv:1503.01808vl [physics.bio-ph] 5 Mar 2015 


Detecting temperature fluctuations at equilibrium 
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Gibbs and Boltzmann definitions of temperature agree only in the macroscopic limit. The ambi¬ 
guity in identifying the equilibrium temperature of a finite sized ‘small’ system exchanging energy 
with a bath is usually understood as a limitation of conventional statistical mechanics. We interpret 
this ambiguity as resulting from a stochastically fluctuating temperature coupled with the phase 
space variables giving rise to a broad temperature distribution. With this ansatz, we develop the 
equilibrium statistics and dynamics of small systems. Numerical evidence using an analytically 
tractable model shows that the effects of temperature fluctuations can be detected in equilibrium 
and dynamical properties of the phase space of the small system. Our theory generalizes statistical 
mechanics to small systems relevant to biophysics and nanotechnology. 


Introduction: Equilibrium properties of a macro¬ 
scopic system exchanging energy with a bath can be de¬ 
scribed by a single intensive paramter, its temperature, 
with remarkable accuracy; independently of the chemi¬ 
cal nature of the bath and system-bath interactions ow¬ 
ing to weak coupling between the system and the bath. 
On the other hand, it is unlikely that a bath couples 
weakly to a system with small number of degrees of free¬ 
dom; consequently, small systems including biophysical 
polymers © and nanomagnets © show considerable de¬ 
viations from the traditional statistical mechanical de¬ 
scription ©. Mathematically, no inverse temperature p 
exists such that the exponential canonical ensemble dis¬ 
tribution accurately predicts equilibrium properties of a 
small system solely dependent on its Hamiltonian. Alter¬ 
natively, Gibbs’ definition of temperature which depends 
on the typical value of energy and the Boltzmann’s defini¬ 
tion of temperature which depends on the mean value of 
energy differ substantially from each other in the case of 
small systems systems 0®. Traditionally, this ambigu¬ 
ity is interpreted as an inevitable statistical uncertainty 
in parameter estimation or a limitation of statistical me¬ 
chanics 0EH9|. 

In this communication, instead of treating the ambigu¬ 
ity in identifying a unique temperature as a limitation, 
we let go of the notion of a unique temperature, espe¬ 
cially for small systems. We identify the ambiguity as 
a consequence of a broad distribution p eq (P)- Further¬ 
more, we identify the broad temperature distribution as 
the f—marginalization of the joint equilibrium distribu¬ 
tion p eq (f, p) of the stochastic variable (f(t),/3(f)) where 
f is the phase space of the system. 

Using maximum entropy arguments, we first estimate 
the joint equilibrium distribution p eq (r,/3) by introduc¬ 
ing two new intensive parameters in the hyperensemble. 
We then show how our theory reduces to traditional sta¬ 
tistical mechanics of macroscopic systems in the suitable 
limit. We illustrate a connections with non-extensive sta¬ 
tistical mechanics of Tsallis mm and our theory at 
thermodynamic equilibrium. Then, we propose Fokker- 
Planck and Langevin equations for the time evolution 
of the instantaneous distribution p(f, P',t). Finally, us¬ 


ing realistic all atom molecular dynamics simulations, we 
present numerical evidence to support our framework and 
discuss its limitations. 

Statitical mechanics of small systems: Consider 
a small system as above. Due to possible non-weak 
coupling between system and the bath, the equilibrium 
phase space distribution of the system p e q (f) will depend 
on the nature of system-bath interactions (fT2HT4| . Let 
us work with the ansatz that the non-canonical behav¬ 
ior arises because the temperature of the system fluc¬ 
tuates © El USD. The joint equilibrium distribution is 
simply 

Peq(r, P) = Peq(f\/3) X p eq {/3). (1) 

In Eq. |T], 

p eq (f\(3) = e^~ H ^ (2) 

is the usual Boltzmann distribution and p e q (/3) needs to 
be determined. Since there are no conservation laws for 
temperature, Gibbs’ ensemble picture is inapplicable. We 
resort to an equally valid alternative. We employ the 
maximum entropy (nraxEnt) framework (16, 11711 . We 
maximize the entropy of the joint distribution p(r, /?) = 
p(r |/3) xp((3) subject to suitable constraints. The entropy 
of the joint distribution is given by 

S \p(r, P)] = ~ X P( f > P) 1o St(g P) ( 3 ) 

f ,/3 

=-^2 p(P) log p(P) + ^2 s(P)p{P) (4) 

0 0 

where 

s(P) = -J2p(r\P) logp(f|/3), (5) 

f 

P(P) = ^2p(r,P) (6) 

f 

is the f— marginal of p(f,P), andp(r|/3) is given by Eq. [2] 

When determining p eq {P) the choice of constraints is 
important. Since the temperature of the system is not 
fixed, we choose (/?) as a constraint. Also, while the en¬ 
tropy of the composite macroscopic system comprising 


2 


the system and the surrounding bath is maximized, the 
entropy of the small system itself not. Consequently, we 
choose the average entropy (s(/J)) as an additional con¬ 
straints and maximize S \p(f, /J)] using Lagrange multi¬ 
pliers. The constraint of average entropy is a common in 
statistical physics and Bayesian statistics of hyperensem¬ 
bles. See OEHZD for different motivations behind this 
choice. After maximization, we find that the equilibrium 
distribution p eq {P) is estimated by 


e As(/3)-C/3 

= -Z{XQ 


(7) 


In Eq. [TJ Z is a generalized partition function and A 
and C, are Lagrange multipliers that determine the shape 
of p e q (/3). If entropy s(/J) is a unitless number, then 
A is unitless and ( has the units of 1//J. The physical 
interpretation of these Lagrange multipliers will become 
clearer below. 

The joint equilibrium distribution p eq (r,p) = 
Peq(r|/?) X Pe q(P) is 


Peq{r,/3) = 


e /3F(/3)-/3ff(r)+As(/3)-C/3 

z [\,0 


( 8 ) 


Thus, instead of describing a thermally equilibrated small 
system with one intensive parameter, its inverse temper¬ 
ature /J, our framework requires two intensive parameters 
A and £ whose meaning will become clear below. 

Connections to traditional statistical mechan¬ 
ics: Assume that the entropy s(/J) is monotonically de¬ 
creasing in /J, a reasonable assumption for systems with 
monotonically increasing density of states. A straightfor¬ 
ward calculation shows that the maximum of p eq (/3) is sit¬ 
uated at /J = Po where /Jo is such that £/A = —c(/Jo)//Jo- 
Here, c(/Jo) is the heat capacity of the system when in¬ 
teracting with an ideal gas at inverse temperature /Jo¬ 
in the limiting case when A —> oo and £ —> oo such 
that their ratio is constant, non-negligible contribution 
to p e q (/J) comes only from near p = po and p eq {P) ~ 
S(/3 — /Jo) where S(x) is the Dirac Delta function. This is 
exactly the traditional canonical ensemble picture where 
the system is assigned the temperature of the surround¬ 
ing thermal bath. It is clear that the magnitudes of A 
and £ dictate the breadth of the p e q (/J) distribution and 
hence the deviation from canonical ensemble. The ratio 
A/C dictates the most likely tempearture of the system. 

Connection to non-extensive statistical me¬ 
chanics: Systems that do not obey the conventional dis¬ 
tributions from statistical mechanics are sometimes en¬ 
tertained within a framework called non-extensive statis¬ 
tical mechanics (Hhd- Though not commonly invoked for 
small systems at equilibrium, here, we will demonstrate 
that non-extensive statistical mechanics can be arrived at 
by marginalization over temperature in a hyperensemble. 

Consider a system whose entropy scales as logarithm 
of temperature, s(/J) = sq log /J, and the internal energy 


scales proportional to the temperature, U(P) = Uo/p, 
when coupled to a bath of ideal gas particles at in¬ 
verse temperature /J. These are excellent assumptions for 
bound systems where density of states increases mono¬ 
tonically with energy. Examples include ideal gas in a 
container and a collection of harmonic oscillators. From 
Eq. [TJ we have 

e _/3 ^/J As °£ Aso+1 


Pe q(/J) — 


r(As 0 +1) 


(9) 


Eq. 22 is a Gamma distribution also known as the gen¬ 
eralized x ~squared distribution. Interestingly, a gamma 
distributed inverse temperature is very commonly used 
in a superstatistical explanation of non-extensive statis¬ 
tics dm. Marginalizing over gamma distributed inverse 
temperature in Eq. [I] results in the so called “Tsallis 
statistics” for the phase space. We have 


Peq{r,P) = 


o P{U(P)-8(P)/P)-f3H{f)+\s{P)-tP 


Z( A,C) 

e Uo-p(H(f)+<;)+(\-l)so log(/3) 


Z(\() 

Integrateing over /J, we have 

p eq (f) oc (1 p 0 (q - 1 )H(f))^ 


( 10 ) 


( 11 ) 


( 12 ) 


Eq. JTTJis the q —generalized canonical ensemble distribu¬ 
tion in Tsallis statistics where 

s o — A A — so +1 , , 

q = ---- and /J 0 = ---. (13) 

So A 1 C 

In the framework of non-extensive statistical mechan¬ 
ics, one arrives at Eq.[l2]by maximizing Tsallis’ q entropy 
with respect to p(r) by constraining an unnatural escort 
expectation of energy m- 

In this work, in contrast to deriving p eq (r) by maximiz¬ 
ing the non-extensive Tsallis entropy by constraining an 
unnatural expectation value, we derive it from a super- 
statistical distribution Eq. [8] and additional assumptions 
about p e q (/J) and system behavior. In our derivation, 
the gamma distribution p eq (/ J) arises in a context specific 
manner i.e. through the logarithmic dependence of the 
entropy on the inverse temperature and by constraining 
average inverse temperature. Therefore, starting from 
the extensive Gibbs-Shannon entropy, maxEnt can act 
as a predictive framework for constructing non-extensive 
effective entropies (El of which the Tsallis entropy is a 
particular example. 

Previously, non-extensive entropies have been criti¬ 
cized from an Occam’s razor point of view dUlMESD 
when compared to the Gibbs-Shannon entropy. Our work 
suggests that non-extensive entropies may arise as ‘ef¬ 
fective entropies’ when considering extensive entropies 
in a hyperensemble. Nevertheless, there is a potential 
loss of information when marginalizing over the temper¬ 
ature P in the hyperensemble that is inherent to con¬ 
structing these effective entropies. We believe that the 
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above demonstration argues in favor the extensive Gibbs- 
Shannon entropy, albeit in a hyperensemble, even when 
the observable phase space may show non-extensive be¬ 
havior. 

Stochastic Dynamics: For simplicity of notation, 
let us consider a one dimensional system. The simplest 
time evolution of the instantaneous distribution p(r, /3; t) 
of the extended phase space that relaxes to a prescribed 
equlibrium distribution p eq (r, /3) can be modeled by an 
over damped Smoluchowski equation. We have 


dp(r,/3; t ) 
dt 


+ D, 


fid Id N 

V7 rdr 7/3 <9/3 y 

d 2 p d 2 p 

.—7 -t-Da—A 


where the ‘forces’ f r and fp are defined as 

d d 

fr = 7^1ogp eq (r,/?) and fi 3 = — logp eq (r, /3). (15) 

By construction, Eq. [14] will relax to the equilibrium dis¬ 
tribution p eq {r, (3) if D r = l/'yr and Dp = l/y p. Note 
that the statistical properties of ( r(t),f3(t )) can also be 
estimated by an overdamped Langevin equation (Brown¬ 
ian dynamics) that is equivalent to Eq. |14| The Langevin 
equation reads 


f = D r f r + \j2D r r\ r 

/3 = Dpfp + y/2Dpr)p (16) 

Here, r/r and rjp are usual uncorrelated Gaussian random 
variables with unit variance. 

Linear analysis: It is instructive to study a linear 
system before analyzing realistic molecules. Consider a 
one dimensional harmonic oscillator interacting with a 
thermal bath. If the deviations from a canonical distri¬ 
bution are negligible, we can treat Eq. [T6] in the linear 
regime by expanding f r and fp to the first order in r and 
/3. In the linear approximation, the joint equilibrium dis¬ 
tribution p eq (r, /3) will be described by a joint normal dis¬ 
tribution. The simplest coupled system of overdamped 
Langevin equations for r(t) and /3(f) that relaxes to to a 
joint normal distribution is given by 

r ~ l\\r + I 121 3 + r] r (17) 

/3 « hir + l 22 /3 + Vp (18) 

We have assumed that the variables r and /3 are appro¬ 
priately scaled by absorbing the diffusion constants D r 
and Dp, lij are the scaled linear expansion coefficients 
of f r and fp, and r] r and r]p are the usual uncorrelated 
Gaussian noises. Integrating over /3(f) and substituting 


in r, we get 

r = lur + h 2 e l22t f ds ■ h 2 ■ e~ l22S 

Jo 

+ l 12 e h2t I ds-r/p ■ e~ l22S + rj r (19) 
Jo 

=> r = (111 + I22 ) r + (^12^21 - I11I22) r 

+ (h 2 V/3 ~ ^22??r-) + Vr (20) 

The time derivative of white noise fj r is a purple noise 
which has quadratically increasing power spectrum. The 
dynamics of temperature fluctuations are governed by the 
linear terms li 2 , l 2 i, l 22 , and the white noise 77 p. These 
terms also appear in the effective Langevin equation for 
r(t). The linear analsysis suggests that one can infer the 
of dynamics of /3(f) by observing the dynamics of r(t). 

The dynamics of r(t) is governed by a much richer 
equation than the usual overdamped Langevin equation. 
A one dimensional small linear harmonic oscillator ex¬ 
changing energy with a thermal bath can be modeled by 
a second order Langevin equation with a combination of 
white and purple noise. These predictions can be tested 
by observing dynamical properties of a small colloidal 
particle trapped in a harmonic well using optical traps. 

A ‘small’ harmonic oscillator: How do we verify 
the effects of temperature fluctuations on the phase space 
of a small system? We resort to realistic molecular dy¬ 
namics simulations of an analytically tractable system 
viz. a harmonic oscillator. 

Consider a three dimensional dumbell shaped Lennard- 
Jones harmonic oscillator interacting non-weakly with a 
bath. Realistic examples include colloidal beads tied to 
each other by a biopolymer or linear molecules such as 
COt. The canonical ensemble distribution for the Har¬ 
monic oscillator is given by 

4/3 3 / 2 r 2 a 2 

Peq (r |/3) = -J^- x e -/3r (21) 

V7T 


where r is the displacement of the oscillator. Without 
loss of generality, we have assumed that the spring con¬ 
stant of the oscillator is k = 2. If the system-bath inter¬ 
actions are non-negligible, we expect that the equilibrium 
phase space distribution of the oscillator will deviate con¬ 
siderably from the Boltzmann distribution. 

The entropy of the oscillator scales as s(/3) ~ log /3 and 
from Eq. [7] we know that the equilibrium distribution 
p e q (/3) will be governed by a Gamma distribution 


Pe q(/3) 


e"^/3 A C A+ i 

r(A + i) 


( 22 ) 


The joint equilibrium distribution p eq (r, /3) = p eq (r\/3) x 
p e q (/3) on the other hand is obtained by multiplying 


Eq. [21] and Eq. |22| 


Peq(r,P) 


4r 2 /3 A+ iC A+1 e _/3 ^ +r2 ) 
V^r(A + i) 


(23) 
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FIG. 1. We study the equilibrium properties of a 3D dumbbell 
shaped harmonic oscillator comprising of Lennard Jones par¬ 
ticles interacting with a bath of water molecules at 300 K (see 
appendix I for details) using all atom MD simulations. The 
numerically obtained marginal distribution p(r) of the oscilla¬ 
tor separation r (black squares) is better captured by Eq. [24 


(red line) than the usual canonical ensemble distribution of 
Eq. 21 (blue line). 


Integrating over all values of /3, we obtain the marginal 
r distribution 


Pe q (r) = 


4r 2 C A+1 r (A + f) (C + r 2 ) 
v^r(A +1) 


-A-J 


(24) 


Moreover, we can also model the dynamics of the oscil¬ 
lator by the coupled Langevin equation of Eq. |16| From 


Eq. 23 the “forces” f r and fp are given by 


fr = ~ ~ 2r/3 and fp = 
r 


3 - 2/3 r 2 - 2/3C + 2A 

w 


(25) 


Eq. 24 along with Eq. 16 where the forces f r and 


are given by Eq. [25] are our predictions for the Harmonic 
oscillator regardless of the bath that is interacting with. 
These predictions can be tested experimentally or in a 
realistic numerical simulation. 

Numerical validation: With the aid of MD simu¬ 
lations of a dumbbell shaped Lennard-Jones harmonic 
oscillator coupled to a bath of water molecules at 300 K 
(see appendix I for details), we confirmed the numerical 
superiority of Eq. [24] compared to Eq. [21] and estimated 
the parameters A ss 2.19 and ( ss 0.34. Fig. |T] shows 
that Eq. [24] which allows for a broad temperature distri¬ 
bution indeed fits the numerically estimated distribution 
much better than the usual canonical ensemble distri¬ 
bution of Eq. [21] It is clear that by allowing the inverse 
temperature to have a broad distribution, the equilibrium 
properties of the harmonic oscillator interacting with its 
thermal surroundings are captured correctly. 


The dynamics of r(t) can be predicted using Eq. 14 by 
studying the equivalent Langevin equation (see appendix 


II). In Fig. [2] we compare the numerically estimated au¬ 
tocorrelation function 


C(r) = (r(r)r(0)) eq - (r)^ q 


(26) 


from MD simulation (black squares) and the prediction 
from the 2-d Langevin equation (red). The predictions 
from an analogous 1-d Langevin equation that relaxes to 
Peq(i") of Eq. 24 are shown in blue. While the dynam¬ 


ics observed in the MD simulation has two time scales 
resulting in a double exponential decay in the autocorre¬ 
lation function, the 1-d Langevin equation is only able to 
capture one effective time scale. On the other hand, the 
2-d Langevin equation has two natural time scales gov¬ 
erned by D r and Dp respectively. The coupled Langevin 
equation equivalent to Eq. [m] (see appendix II) with 
dt = 5 x 10” 8 and Dp ss 50 x D r does indeed captures the 
autocorrelation function while an analogous 1-d equation 
fails to do so (see appendix II for details of the fit). 

In appendix III we show that the theoretical predic¬ 
tions are valid over a range of bath temperatures and 
system-bath interactions. In this work, we study a sys¬ 
tem whose canonical ensemble distribution can be analyt¬ 
ically computed and the entropy analytically estimated. 
This allowed us to compute p e q(/3) and p eq (r,/3) analyt¬ 
ically. For more realistic systems with multiple degrees 
of freedom, p eq (r\/3) needs to be estimated numerically 
along with p eq {/3)- 

In summary, the mesoscopic harmonic oscillator inter¬ 
acting with a thermal bath of water molecules shows sig¬ 
nificant deviation from the canonical ensemble descrip¬ 
tion. We can correctly predict both equilibrium and dy¬ 
namical properties of the oscillator by allowing its tem¬ 
perature to vary as a stochastic variable which is coupled 
with the phase space variable r(t). 


Discussion: 

It is known that, at equilibrium, mesoscopic systems 
have larger fluctuations compared to a macrosopic sys¬ 
tem. We have argued that these enhanced fluctuations 
be understood as arising from a dynamically fluctuating 
temperature. 

How do we reconcile a time dependent temperature, 
a non-equilibrium phenomena prima facie, in an equilib¬ 
rium setting? Even though the temperature is chang¬ 
ing, the extended phase space (r(t),/3(t)) is still gov¬ 
erned by a detailed-balanced Markov process. It’s an 
easy calculation to show that the entropy production, 
as defined in stochastic thermodynamics (E3, is indeed 
zero for the hyperensemble. Nevertheless, there are mul¬ 
tiple questions which need resolution. For example, How 
do we formulate non-equilibrium phenomena in the hy¬ 
perensemble setting? For example, how do we modify 
non-equilibrium fluctuation relationships (128|) for small 
systems? We leave this to future work. 
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FIG. 2. We study the autocorrelation function of the har¬ 
monic oscillator interacting with a bath of water molecules. 
We model the dynamics of the extended phase space 
(r(t), /3(t)) using a simple coupled Langevin equation (see ap¬ 
pendix II). We find that the 2 dimensional Langevin equation 
(red line) captures the two time scales inherent to the dynam¬ 
ics of r(t) as observed in MD simulations (black squares). On 
the other hand, an analogous 1-d Langevin equation can only 
capture one effective time scale (blue line). 
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APPENDIX I: MD SIMULATIONS 


A harmonic dumbbell oscillator consisting of two 
Lennard-Jones particles was immersed in a bath of 333 
TIP3 Pill I5U[) water molecules. NVT molecular dynam¬ 
ics simulations were run with NAMD pij) at 300K with 
a box size of 19.12A. The CHARMM (l32|) forcefield was 
used to describe the interaction between the harmonic os¬ 
cillator particles and surrounding water molecules. The 
spring constant for the dumbell was chosen to be k = 
0.25 kcal/mol-A 2 , the e parameter was set at e = —20.0 
kcal/mole and the size parameter was set at r = lA. 
The systems were minimized for 2000 steps followed by 
an equilibration of 1 nanosecond and a production run of 
2 nanosecond. The integration time step was 0.25 fem¬ 
toseconds and the trajectory was saved every 2.5 fem¬ 
toseconds. 


APPENDIX II: FITTING LANGEVIN 
DYNAMICS TO DATA 


The coupled Langevin equation corresponding to 
where the equilibrium distribution p eq (r, /?) is 


14 


Eq. _ 
given by Eq. [22] is given by 


/ r(t + dt) 
P(t + dt) 


( r{t) 



+ V2dt 


T] r 's/~Dr A 
r]fty/D~0 ) 


(27) 


Here, r] r and rjg are uncorrelated Gaussian random vari¬ 
ables with unit variance, dt is a small time step, D r and 
Dp are diffusion coefficients for the phase space coordi¬ 
nate r and the temperature /3. 

From the MD simulation, we first estimated the auto- 
correlaion function C(t). The Langevin equation can be 
scaled in time by multiplying the diffusion constants and 
dividing the time step dt by the same number. In order 
to ensure smooth integration, we first set the integration 
time step to a very small value; dt = 5 x 10 -8 . Every 
pair ( D ri Dp ) of diffusion constants predicted an auto¬ 
correlation function that had two inherent time scales 
manifested in a double exponential decay. We man¬ 
ually scanned the (D r ,Dp )—space to match the MD- 
autocorrelation function. We found that D r = 1 and 
Dp = 50 gave reasonable fits (red curve). 

We also wrote down a 1-d Langevin equation analogous 
to Eq. [27j 


r(t + dt) « r(t) + D r f r dt + \j2D r dtr\ r 


(28) 


where f r = j^\ogp eq (r) (see Eq. 24). This equation 
had only one diffusion constant D r . A one dimensional 
scan of D r suggested that the autocorrelation function 
predicted using the 1-d Langevin equation always had a 
single exponential decay. We found the best fit to the 
autocorrelation function at D r ss 50 (blue curve). 





